##script adapted from http://www.adamlilith.net/!ecology/sdmShortCourseKState2012/grinnellExercise/exercise-tuning-maxent-using-beta-and-aic.html
##steps 7-10 are the same from script: Exercise III - Tuning Maxent Using Beta and AIC.r

source('/home/jc148322/scripts/libraries/cool_functions.r')
library(SDMTools)

timeline='Pre_decline' # or 'Pre_decline'

wd= paste('/home/jc148322/frog_declines/SDM_AICC/',timeline,'/models/',sep='') #your /models/ directory

species=list.files(wd)
out=NULL


for (spp in species) {
	spp.dir = paste(wd,spp,'/',sep=''); setwd(spp.dir) #set the working directory to the species directory
		
	var.dirs=list.files(spp.dir)
	var.dirs=grep('bioclim',var.dirs,value=TRUE)
	var.dirs=var.dirs[rev(order(sapply(strsplit(var.dirs,'\\.'),length)))] #reorder by decreasing model complexity
	
	AICCs=AICs=kcount=NULL
	#cycle through each of the species
	for (v in var.dirs) { cat(v,'\n')
		var.dir=paste(spp.dir,v,'/',sep='');setwd(var.dir)
		occur=read.csv('occur.csv',as.is=TRUE)
		#read in ascii of raw values
		tasc=read.asc.gz('output/ascii/current.76to05.asc.gz')
		
		### STEP 7 extract predictions at species' presence sites
		values=extract.data(cbind(occur$lon,occur$lat),tasc)
		
		### STEP 8 calculate log likelihood
		L <- sum( log(values), na.rm=TRUE) # calculate log likelihood

		### STEP 9 calculate number of parameters used by Maxent model
		# The number of parameters is obtained from the "lambdas" file
		Lambdas=readLines(paste('output/',spp,'.lambdas',sep=''))
		
		k <- 0 # number of parameters

		for (thisLambda in Lambdas <- Lambdas) { # for each line in lambda object

			commaPos <- gregexpr( text=thisLambda, pattern=',') # get location of commas
			
			if (length(commaPos[[1]]) > 1) { # if there is >1 comma in this line (this is not a parameter line)
				paramValue <- as.numeric( substr(x=thisLambda, start=commaPos[[1]][1]+1, stop=commaPos[[1]][2]-1 ) ) # convert string between first two commas to numeric
				if (paramValue !=0) k <- k + 1 # increment number of parameters
			} # if there is >1 comma in this line
		}

		## STEP 10 calculate AICc
		aicc <- -2 * L + 2 * k + (2 * k * (k + 1)) / ( nrow(occur) - k - 1) # AICc
		aic = -2 * L + 2 * k
		
		AICCs=c(AICCs, aicc)
		AICs=c(AICs, aic)
		kcount=c(kcount,k)
	}
	out=cbind(var.dirs,round(AICCs,2),AICs, kcount)
	colnames(out)=c('variables','AICc values','AIC values','#parameters')
	
	setwd(spp.dir)
	write.csv(out,'aicc_values.csv',row.names=FALSE)


}

# The model with the lower value of AICc is teh most parsimonious, so its beta should be used in subsequent moding of this species, assuming the training records and predictors do not change.  These steps can be perfomed using other values of beta (be sure to change variable names if you do!) and could be done in a loop over multiple values of beta.
